Reconstructing the Density of States by History-Dependent Metadynamics 
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We present a novel method for the calculation of the energy density of states D(E) for systems described by 
classical statistical mechanics. The method builds on an extension of a recently proposed strategy that allows 
the free energy profile of a canonical system to be recovered within a pre-assigned accuracy, [A. Laio and M. 
Parrinello, PNAS 2002]. The method allows a good control over the error on the recovered system entropy. 
This fact is exploited to obtain D(E) more efficiently by combining measurements at different temperatures. 
The accuracy and efficiency of the method are tested for the two-dimensional Ising model (up to size 50x50) by 
comparison with both exact results and previous studies. This method is a general one and should be applicable 
to more realistic model systems. 
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It has long been recognized that all energy -related ther- 
modynamic properties of a classical canonical system can 
be calculated once the energy density of states D(E) is 
known. In fact, starting from the partition function Z = 
J dE D(E)e~ /3E ', quantities like free energies, specific heats 
and phase transition temperatures can be computed in a 
straightforward manner. In principle, all canonical aver- 
ages can be calculated through a multidimensional density of 
states. Due to this central role in equilibrium thermodynam- 
ics a variety of theoretical and computations studies have ad- 
dressed the problem of how to obtain D(E) (or, equivalently, 
the entropy profile S(E) — \n(D(E)) in a reliable and effi- 
cient way. 

In principle, D(E) could be calculated from the histogram 
of the energies visited during a single "very long" constant- 
temperature simulationl 1]. In practice, for any finite simu- 
lation only a limited energy windows is sampled so that the 
recovery of the system thermodynamics over a wide tempera- 
ture range is unfeasible. 

Several alternative strategies have been developed to rem- 
edy this shortcoming. For example the multiple histogram 
reweighting technique relies on performing several simula- 
tions at different temperatures, so as to explore different (over- 
lapping) energy intervals |2]. The various histograms are then 
optimally combined to obtain D(E) over the union of the en- 
ergy intervals. Another successful family of techniques aims 
at obtaining D(E) by changing iteratively the probability with 
which the various energy levels are visited in stochastic dy- 
namics until the recorded energy histogram is "flat". Such 
methods include entropic sampling JiLmulticanonical and 
broad-histogram techniques techniques |3, 5] and also the re- 
cent method of Wang and Landau ||6[ . 

These and similar techniques have proved to be very valu- 
able in a variety of contexts |7, 8, 9], but there is still ample 
scope for alternative approaches that could provide improved 
efficiency and better error control. Here we propose to mod- 
ify and extend the metadynamics method recently introduced 
by two of us 11011 for evaluating D(E). The algorithm we 
introduce allows a good error control on the explored energy 
range. Moreover, although within our approach D(E) could 
be reconstructed performing simulations at virtually any tem- 



perature, it is still possible to exploit the Boltzmann bias to 
focalize the computational effort for exploring the region of 
phase space of relevance for a temperature of interest (e.g. of 
a phase transition). In the following we will first describe this 
general strategy followed by its application to the typical ref- 
erence case constituted by the two-dimensional Ising model. 

The essence of the metadynamics approach is first to iden- 
tify those relevant collective variables (CVs) that are difficult 
to sample. In other applications lHoll a similar metodology 
was applied in order to observe a specific transition (e.g. a 
chemical reaction) in systems described by a complex atom- 
istic Hamiltonian. At this scope CVs explicitly depending on 
the microscopic configuration of the system have been em- 
ployed, like, e.g., coordination numbers or distances between 
specific atoms of the system. Here, we aim at reconstructing 
the canonical free-energy profile, F(E) = E — T S(E) and 
the relevant variable is E (which is also an explicit function of 
the microscopic configuration of the system). At each meta- 
dynamics step the system evolution is guided by the gener- 
alised force which combines the action of the thermodynamic 
force (which would trap the system in free-energy minima) 
dF(E)/dE and a history-dependent one which disfavours 
system configurations already visited. The history-dependent 
potential, Fq, is constructed as a sum of Gaussians of width 
AE and height w centred around each value of E already ex- 
plored during the dynamics. As shown in ref 1 10], Fq fills in 
time the minima in the free energy surface and, in the limit of 
a long metadynamics F(E) + Fq{E) tends to become flat as 
a function of E and hence —Fg(E) becomes an approximant 
of F(E). 

Clearly, the exact form of the metadynamics equations and 
the choice of the parameters w and AE may affect signifi- 
cantly the accuracy of this estimate. Furthermore careful con- 
trol of the error is essential for reconstructing a reliable den- 
sity of states. This requires that the algorithm in ref. lHoll 
be substantially improved. The modified metadynamics equa- 
tions are: 
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where \ is a random number between and 1 . The general- 
ized force f(E) is given by / (E) = + 
with the history-dependent potential, i 7 ^, defined as 

Fq {E)=y2,we~ i aA J } . (3) 

tt<t 

By displacing the center of the Gaussian -E^ 1 " 1 with respect to 
the point of evaluation of the generalised force (equation^, 
the added Gaussian maximally compensates f{E t ) and flat- 
tens F + Fq around E*. The energy step performed at ev- 
ery metadynamics iteration is chosen randomly in the interval 
between AE and 1.5Ai? (equation|2j in order to reduce the 
correlation induced by the dynamics. 




FIG. 1: Example of reconstruction of a preassigned (analytic) free- 
energy profile (black curve), F(E), by means of metadynamics 
runs consisting of 200 Gaussians of spread AE = 0.4 and height 
w = 0.16. To mimic the uncertainty on the thermodynamic force 
encountered in stochastic simulations, a Gaussian noise of width 0.3 
was added to F'(E). The broken red and indigo lines denote the 
filled profile, SF(E) = F(E) - Fr(E), obtained respectively with 
an unsmoothed (r c = 0) and smoothed (r c = 100) metadynam- 
ics. Notice the shorter correlation lengths and smaller amplitudes 
of the fluctuations in the second case. The average and dispersion 
of 5F(E) were calculated (blue and brown lines wth error bars) 
for 1000 independent runs, revealing both the absence of biases in 
Fr(E) and the constancy of the error except close to the boundaries. 
By changing w and modifying the number of Gaussians so as to work 
at a constant filled free-energy volume, we calculated over the range 
\E\ < 2.5 the average, maximum and minimum values of SF(E) 
measured over 1000 runs, see inset. The linear dependence of the 
error on Fr is thus apparent both for the unsmoothed and smoothed 
metadynamics (red and blue respectively). 

Finally, in order to further reduce the spatial correlations 
in Fq, we notice that when the metadynamics is terminated, 



say at position E l , Fq will present a bump in a region around 
E l whose spread depends on the correlation time of the meta- 
dynamics. This effect can be lessened if the contribution of 
the Gaussians placed at the end of the dynamics are weighted 
less. Therefore, after the metadynamics with constant w is 
terminated we reconstruct the free energy from Fr (E) = 

— J2 u <t w ' an ' 1 ( ) e 2&Ei where r c is taken to be 
larger than the typical time required to sweep the "filled" en- 
ergy range. Other smoothing functions, such as min(l, 
can, of course, also be chosen and have an analogous effect. 

The modified metadynamics algorithm allows the efficient 
reconstruction of F(E) in the explored energy range, within 
an error that is ultimately controlled only by the Gaussian's 
height w. We demonstrate the quality of the algorithm by re- 
constructing an "ideal", pre-assigned, F (E) (see Fig. [Q. If 
the method is void of systematic biases one would expect the 
quantity SF (E) = Fr (E) — F (E) to be, on average, con- 
stant throughout the filled energy range. Moreover, we would 
expect deviations from the constant average value to be of the 
order of w. These properties are confirmed by the results pre- 
sented in the inset of Fig. 1, where the uniformity of the aver- 
age value of 8F(E) is apparent, together with the constancy 
of its dispersion, a 2 (E) = {SF (E) 2 ) - (6F (E)) 2 . The plots 
also illustrate the benefits of the "smoothing" procedure over 
the last part of the metadynamics trajectory, since this results 
in a decrease of the spread of 5F(E). As is visible in the inset 
of Fig. 1, the dispersion is further confirmed to be approxi- 
mately proportional to w. An important fact is that near the 
boundaries of the explored energy interval Fr(E) decays to 
zero and hence deviates from the true free energy. To identify 
the interval over which Fr is reliable we need to ascertain if 
the number of Gaussians accumulated at a given energy E, 
« \Fr(E)/w\, is significantly larger than the minimal num- 
ber of superposed Gaussians needed to produce the observed 
free energy derivative, w F 1 (E)/(w/AE). In this work we 
have required that the ratio of the former to the latter be greater 
than 5. 

We now use the algorithm described above to compute 
D(E) for an NxN two-dimensional Ising model with ferro- 
magnetic nearest-neighbour interactions and periodic bound- 
ary conditions In fact, exact expressions for S(E) are 
available 1 12] and the error induced by the algorithm can be 
explicitly estimated and comparison with other approaches 
can be made[6|. 

Virtually the entire computational effort of the metadynam- 
ics is spent in estimating the thermodynamic forces dF/ dE 
at each energy value, E, visited by the metadynamics. To 
respect the discrete nature of the system's energy spectrum, 
the continuous value of E 1 produced by eq. (0 was dis- 
cretized to the nearest energy level. Due to the discreteness 
of the energy spectrum, the thermodynamic force in E cannot 
be estimated by the Lagrange-multipliers technique described 
in Ref.| 10], and is rather obtained using a centred difference 
approach. If pi and p2 are the occupation probabilities of 
the two energy levels, E% and E2 adjacent to E (we assume 
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\E - E-A = \E- E 2 \), we have: 
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pi and p2 are evaluated with an umbrella-sampling strategy 
consisting of a Monte Carlo evolution of the system (Metropo- 
lis acceptance/rejection of single-spin flips) under the action 
of an effective Hamiltonian obtained by adding to the energy 
of a given spin configuration, E, the term 1/2 K(E — E) 2 . 
A suitable choice of the parameter K forces the system to 
explore the energy region around E thus populating appre- 
ciably both levels E\ and E 2 . The symmetry with respect 
to E of the added umbrella potential allows to calculate the 
thermodynamic force through the same eq. |4] with pi and 
P2 being the fraction of times that the corresponding energy 
levels are encountered in n statistically independent configu- 
rations picked with the modified canonical weight. For this 
purpose the Monte Carlo trajectory was sampled at intervals 
comparable to the autocorrelation time after having discarded 
a few tens of initial system sweeps. The MC sampling was 
stopped when the estimated uncertainty on the force II 311 was 
equal to the maximum force introduced by a single Gaussian, 
w exp(— 1/2) /AE. This choice ensures that, for large values 
of t, f(E) is of the order of w/AE. If the force is calculated 
with much greater accuracy the repeated superposition of the 
Gaussians would still lead to an uncertainty of order w/AE 
on f(E). 

By means of such a metadynamics it is therefore possi- 
ble to reconstruct the free energy profile, Fr(E, T). An es- 
timator for the system entropy is given by Sr(E) = [E — 
F R (E, T)]/T. The uncertainty over F R (E, T) is inherited by 
Sr(E) whose a priori dispersion is thus of the order w/T. 
Thus, the expected error on the entropy profile is constant. 
This represents a major difference over standard reweighting 
techniques, where the accuracy on the calculated entropy usu- 
ally deteriorates as one moves away from the free energy min- 
ima. 

If the goal is to reconstruct S(E) over a wide range of en- 
ergy, it seems natural to combine the outcome of several meta- 
dynamics at various temperatures, in analogy with multiple- 
histogram techniques|2|. In the following we shall indicate 
with Srj the entropy reconstructed in the ith metadynamics 
carried out at temperature Ti and with Gaussians of height 
Wi and width 5i. Due to the temperature-dependence of the 
free energy, each run will typically explore a different energy 
range. The data obtained in the different metadynamics runs 
can be optimally combined to provide a single entropy esti- 
mate, S, over the union of the explored regionsQ. To do so 
we recall that the entropy Sn.i is known only up to an addi- 
tive constant Cj and that the uncertainty on Srj is ti — Wi/Ti 
throughout the reliably-explored energy range. 
This leads us naturally to consider a maximum likelihood ap- 
proach to obtain S and the additive constants by minimizing 
the least-squares function, 
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FIG. 2: Results from the metadynamics runs carried out on the 32x32 
Ising system. Top: the exact and reconstructed entropies. In both 
cases the entropy has been normalised so that ^2 E exp[S(E)] = 1. 
The abscissa indicates the energy per spin. Middle: Difference be- 
tween the true and reconstructed entropy. The dispersion expected a 
priori on the reconstructed entropy is shown with a red line. Bottom: 
the exact and reconstructed specific heat and average energy (inset) 
as a function of temperature. 
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where the first sum is carried over the various metadynamics 
runs and the second one runs over the system (discrete) en- 
ergy levels with the proviso that ti is equal to infinity outside 
the reliable energy range. The determination of S(E) through 
the minimization of C relies on the statistical independence of 
each term in the sum of equation This is realized only ap- 
proximately due to the existence of an intrinsic scale of auto- 
correlation for the reconstructed free-energy/entropy dictated 
by the correlation lenght of the metadynamics. Therefore, the 
minimization of is meaningful provided that each energy 
value is covered by several metadynamics runs, each explor- 
ing an interval substantially larger than the Gaussian widths. 
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The requirement of stationarity for C leads to self- 
consistent equations which can be solved iteratively in terms 
of the Cj's and S(E). Despite the presence of the ad- 
ditive terms, Cj's, which distinguish the present juxtaposi- 
tion scheme from others already available, the self-consistent 
equations are simple both in their formulation and numeri- 
cal implementation. Convergence to the solutions is typically 
reached in a few tens of iterations. The least-squares approach 
also allows the expected standard deviation of S(E) to be cal- 
culated: 

a(S(E))-^ = J2 ^ E r 2 ■ ( 6 ) 

i—l..n 

The use of the reconstruction strategy is first illustrated for 
the 32x32 Ising model, see Fig. The curve for the entropy, 
Fig. |21 results from the combination of runs at six tempera- 
tures, T — 2, 2.6, 3.0, 3.4, 6.0 and 12.0. At each temperature 
1000 Gaussians where used with w/T w 0.5, AE w 0.04, 
K w 0.4 and r c = 300. The total computational effort re- 
quired 7.5 10 5 MC sweeps. The comparison with the true 
system entropy reveals that S(E) was correctly reconstructed 
throughout the explored energy range [-1.93; 1.93] (exploit- 
ing the ferro/antiferro symmetry) with an uncertainty that is 
approximately constant (its average being 0.17) and in agree- 
ment with the one expected a priori from eq. |6] The corre- 
sponding average relative error on S(E) was 0.05%, similar 
to that obtained with a comparable number of sweeps in the 
recent and powerful approaches described in refs. ||5j,|6j|. 

Analogous runs were repeated for the 50x50 system using 
three temperatures, T — 2.6,4 and 12 and the same param- 
eters as before. By using 2.2 10 6 Monte Carlo sweeps S(E) 
was reconstructed over the energy range [-1.8,1.8] again with 
an average error of 0.24, again in agreement with the one ex- 
pected a priori. This confirms that the proposed strategy al- 
lows good control over the final accuracy on S. 

We wish to point out that, in the high-temperature limit, 
our approach has strong analogies with the Wang and Lan- 
dau algorithm, in which D(E) is also modified in a history- 
dependent fashion|6] (their "pointwise" modification of the 
density of states can be viewed as a limiting case of our Gaus- 
sians). As in their case, with one metadynamics run at a single 
temperature we could explore the whole energy range. This, 
however, may be inefficient, especially in a realistic model 
since it could require an impractically large number of Gaus- 
sians. Within our approach, it is not necessary to renounce to 
the Boltzmann bias, and it is possible to focalize the effort for 
exploring the region of phase space of relevance for the tem- 
perature of interest. With this respect, our metodology can 
be viewed as a finite temperature extension of the Wang and 
Landau algorithm. 

Although in its present formulation the proposed method 



allows an accurate and efficient recovery of a system entropy, 
it is certainly susceptible to further generalizations and im- 
provements. In particular, in order to improve the resolution, 
the height and/or width of the Gaussians may be changed as 
the metadynamics progresses, in analogy with the method of 
ref. |0. The application of the method to first order phase 
transitions is conceptually straightforward although, in prac- 
tice, the elimination of hysteretic effects in the metadynam- 
ics may prove computationally expensive. However, these ef- 
fects could be eliminated by exploiting the ability of metady- 
namics to sample multidimensional free energy surfaces lllOII . 
Supplementing E with auxiliary order parameters suitable for 
characterizing the transition should facilitate the overcoming 
of the free energy barriers associated with the nucleation of 
the new phase and thus eliminate/reduce the hysteresis. The 
progress made here constitutes a substantial improvement to 
the accuracy of the metadynamics approach and illustrates 
its relation to other very powerful methods like multiple his- 
togram reweightingfl and Wang and Landau algorithm|6|. 
Given the potential range of applications of metadynamics we 
expect that our work will have an impact far broader than the 
present demonstrative calculation on the Ising model. 
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